home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpkdrv.zip / SUD.FOR < prev    next >
Text File  |  1984-04-30  |  13KB  |  440 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)  This is doesn't work on PC's!
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SQRXX(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SQRXX(LUNIT)
  14.       INTEGER LUNIT
  15.       REAL RX(20,10),R(10,10),XROW(10),Z(10,4),YROW(10),S(10)
  16.       REAL SCALE,TINY,HUGE,RHO(2),C(10),SMACH
  17.       INTEGER N,P,LDX,LDR,LDZ,CASE
  18.       LOGICAL NOTWRT,OFLOW,UFLOW
  19.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  20.       LDZ = 10
  21.       LDX = 20
  22.       LDR = 10
  23.       NOTWRT = .TRUE.
  24.       OFLOW = .FALSE.
  25.       UFLOW = .FALSE.
  26.       TINY = SMACH(2)
  27.       HUGE = SMACH(3)
  28.       SCALE = 1.0E0
  29.       DO 60 CASE = 1, 3
  30.          GO TO (10, 20, 30), CASE
  31.    10    CONTINUE
  32.             N = 20
  33.             P = 10
  34.          GO TO 40
  35.    20    CONTINUE
  36.             N = 10
  37.             P = 4
  38.          GO TO 40
  39.    30    CONTINUE
  40.             N = 10
  41.             P = 1
  42.    40    CONTINUE
  43.          CALL SGETRX(RX,LDX,N,P,S)
  44.          WRITE (LUNIT,130) CASE,N,P
  45.          IF (NOTWRT) GO TO 50
  46.             WRITE (LUNIT,160)
  47.             CALL SARRAY(RX,LDX,N,P,P,LUNIT)
  48.    50    CONTINUE
  49.          CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  50.          CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  51.          CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  52.    60 CONTINUE
  53.       CASE = 4
  54.       N = 10
  55.       P = 4
  56.       WRITE (LUNIT,130) CASE,N,P
  57.       WRITE (LUNIT,140)
  58.       CALL SGETRX(RX,LDX,N,P,S)
  59.       DO 80 J = 1, P
  60.          DO 70 I = 1, N
  61.             RX(I,J) = HUGE*RX(I,J)
  62.    70    CONTINUE
  63.    80 CONTINUE
  64.       SCALE = HUGE
  65.       OFLOW = .TRUE.
  66.       IF (NOTWRT) GO TO 90
  67.          WRITE (LUNIT,160)
  68.          CALL SARRAY(RX,LDX,N,P,P,LUNIT)
  69.    90 CONTINUE
  70.       CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  71.       CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  72.       CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  73.       OFLOW = .FALSE.
  74.       CASE = 5
  75.       N = 10
  76.       P = 4
  77.       WRITE (LUNIT,130) CASE,N,P
  78.       WRITE (LUNIT,150)
  79.       CALL SGETRX(RX,LDX,N,P,S)
  80.       DO 110 J = 1, P
  81.          DO 100 I = 1, N
  82.             RX(I,J) = TINY*RX(I,J)
  83.   100    CONTINUE
  84.   110 CONTINUE
  85.       SCALE = TINY
  86.       UFLOW = .TRUE.
  87.       IF (NOTWRT) GO TO 120
  88.          WRITE (LUNIT,160)
  89.          CALL SARRAY(RX,LDX,N,P,P,LUNIT)
  90.   120 CONTINUE
  91.       CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  92.       CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  93.       CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  94.       RETURN
  95. C
  96.   130 FORMAT (11H1    CASE =, I2, 5X, 3HN =, I2, 5X, 3HP =, I2 /////)
  97.   140 FORMAT (22H         OVERFLOW TEST /////)
  98.   150 FORMAT (24H          UNDERFLOW TEST /////)
  99.   160 FORMAT (5H   RX)
  100.       END
  101.       SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT)
  102.       INTEGER LDA,M,N,NNL,LUNIT
  103.       REAL A(LDA,1)
  104. C     X
  105. C          FORTRAN IABS,MIN0
  106. C
  107.       INTEGER I,J,K,KU,NL
  108.       NL = IABS(NNL)
  109.       IF (NNL .LT. 0) GO TO 30
  110.          DO 20 I = 1, M
  111.             WRITE (LUNIT,70)
  112.             DO 10 K = 1, N, NL
  113.                KU = MIN0(K+NL-1,N)
  114.                WRITE (LUNIT,70) (A(I,J), J = K, KU)
  115.    10       CONTINUE
  116.    20    CONTINUE
  117.       GO TO 60
  118.    30 CONTINUE
  119.          DO 50 J = 1, N
  120.             WRITE (LUNIT,70)
  121.             DO 40 K = 1, M, NL
  122.                KU = MIN0(K+NL-1,M)
  123.                WRITE (LUNIT,70) (A(I,J), I = K, KU)
  124.    40       CONTINUE
  125.    50    CONTINUE
  126.    60 CONTINUE
  127.       RETURN
  128. C
  129.    70 FORMAT (1H , 4(2E13.6, 4X))
  130.       END
  131.       SUBROUTINE SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  132.       INTEGER LDR,LDX,LDZ,P,LUNIT
  133.       REAL R(LDR,1),RX(LDX,1),Z(LDZ,1),XROW(1),YROW(1),S(1)
  134.       REAL RHO(1),SCALE,C(1)
  135.       LOGICAL NOTWRT,OFLOW,UFLOW
  136.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  137.       INTEGER I,J,JP1,PM1
  138.       WRITE (LUNIT,50)
  139.       CALL SCHDD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S,INFO)
  140.       PM1 = P - 1
  141.       IF (PM1 .LT. 1) GO TO 30
  142.       DO 20 J = 1, PM1
  143.          JP1 = J + 1
  144.          DO 10 I = JP1, P
  145.             R(I,J) = 0.0E0
  146.    10    CONTINUE
  147.    20 CONTINUE
  148.    30 CONTINUE
  149.       IF (NOTWRT) GO TO 40
  150.          WRITE (LUNIT,80)
  151.          CALL SARRAY(R,LDR,P,P,P,LUNIT)
  152.          WRITE (LUNIT,60)
  153.          CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
  154.          WRITE (LUNIT,70) (RHO(I), I = 1, 2)
  155.    40 CONTINUE
  156.       CALL SMPDD(R,LDR,P,RX,LDX,Z,Z(1,3),LDZ,RHO,LUNIT)
  157.       RETURN
  158.    50 FORMAT ( /////
  159.      *         46H     STEP THREE    DOWNDATING XROW,YROW AND Z, ///)
  160.    60 FORMAT ( /// 4H   Z)
  161.    70 FORMAT ( /// 6H   RHO // 1H , 2E13.6)
  162.    80 FORMAT (4H   R)
  163.       END
  164.       SUBROUTINE SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  165.       INTEGER N,P,LDR,LDX,LUNIT
  166.       REAL R(LDR,1),RX(LDX,1),XROW(1),S(1)
  167.       REAL C(1)
  168.       LOGICAL NOTWRT,OFLOW
  169.       REAL RELM,XELM,Y,Z
  170.       REAL XMAX,RMAX,TEST,T,SCALE,SMACH
  171.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  172.       DO 20 I = 1, P
  173.          DO 10 J = 1, P
  174.             R(I,J) = 0.0E0
  175.    10    CONTINUE
  176.    20 CONTINUE
  177.       DO 40 I = 1, N
  178.          DO 30 J = 1, P
  179.             XROW(J) = RX(I,J)
  180.    30    CONTINUE
  181.          CALL SCHUD(R,LDR,P,XROW,Z,1,0,Y,Y,C,S)
  182.    40 CONTINUE
  183.       WRITE (LUNIT,100)
  184.       IF (NOTWRT) GO TO 50
  185.          WRITE (LUNIT,120)
  186.          CALL SARRAY(R,LDR,P,P,P,LUNIT)
  187.    50 CONTINUE
  188.       RMAX = 0.0E0
  189.       XMAX = 0.0E0
  190.       DO 90 I = 1, P
  191.          DO 80 J = 1, I
  192.             RELM = 0.0E0
  193.             XELM = 0.0E0
  194.             NIM = MIN0(I,J)
  195.             DO 60 K = 1, NIM
  196.                RELM = RELM + R(K,I)/SCALE*(R(K,J)/SCALE)
  197.    60       CONTINUE
  198.             DO 70 K = 1, N
  199.                XELM = XELM + RX(K,I)/SCALE*(RX(K,J)/SCALE)
  200.    70       CONTINUE
  201.             T = AMAX1(ABS(XELM),0.0E0)
  202.             XMAX = AMAX1(XMAX,T)
  203.             T = AMAX1(ABS(XELM-RELM),0.0E0)
  204.             RMAX = AMAX1(RMAX,T)
  205.    80    CONTINUE
  206.    90 CONTINUE
  207.       TEST = RMAX/XMAX/SMACH(1)
  208.       WRITE (LUNIT,110) TEST
  209.       RETURN
  210. C
  211.   100 FORMAT ( ///// 25H    STEP ONE   UPDATING X ///)
  212.   110 FORMAT ( /// 15H     STATISTICS //
  213.      *         42H      RH*R    ............................, E10.2)
  214.   120 FORMAT (4H   R)
  215.       END
  216.       SUBROUTINE SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  217.       INTEGER LDR,LDX,LDZ,P,LUNIT
  218.       REAL R(LDR,1),RX(LDX,1),XROW(1),YROW(1),Z(LDZ,1),S(1)
  219.       REAL RHO(1),C(1)
  220.       REAL SCALE
  221.       LOGICAL NOTWRT,OFLOW,UFLOW
  222.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  223.       DO 20 I = 1, P
  224.          DO 10 J = 1, P
  225.             RX(I,J) = R(I,J)
  226.    10    CONTINUE
  227.    20 CONTINUE
  228.       DO 40 I = 1, P
  229.          Z(I,1) = 0.0E0
  230.          DO 30 J = I, P
  231.             Z(I,1) = Z(I,1) + R(I,J)
  232.    30    CONTINUE
  233.          Z(I,2) = Z(I,1) + SCALE
  234.          Z(I,3) = Z(I,1)
  235.          Z(I,4) = Z(I,2)
  236.    40 CONTINUE
  237.       DO 60 I = 1, P
  238.          XROW(I) = 0.0E0
  239.          DO 50 J = 1, I
  240.             XROW(I) = XROW(I) + R(J,I)
  241.    50    CONTINUE
  242.    60 CONTINUE
  243.       YROW(1) = 0.0E0
  244.       DO 70 I = 1, P
  245.          YROW(1) = YROW(1) + XROW(I)
  246.    70 CONTINUE
  247.       YROW(2) = YROW(1) - SCALE
  248.       RHO(1) = SCALE
  249.       RHO(2) = SCALE
  250.       IF (NOTWRT) GO TO 80
  251.          WRITE (LUNIT,100)
  252.          CALL SARRAY(XROW,P,P,1,-P,LUNIT)
  253.          WRITE (LUNIT,110)
  254.          CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
  255.          WRITE (LUNIT,120)
  256.          CALL SARRAY(YROW,2,2,1,-2,LUNIT)
  257.    80 CONTINUE
  258.       CALL SCHUD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S)
  259.       WRITE (LUNIT,130)
  260.       IF (NOTWRT) GO TO 90
  261.          WRITE (LUNIT,150)
  262.          CALL SARRAY(R,LDR,P,P,P,LUNIT)
  263.          WRITE (LUNIT,110)
  264.          CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
  265.          WRITE (LUNIT,140) (RHO(I), I = 1, 2)
  266.    90 CONTINUE
  267.       CALL SMPUD(R,LDR,RX,LDX,P,XROW,Z,LDZ,RHO,LUNIT)
  268.       RETURN
  269. C
  270.   100 FORMAT ( ///// 7H   XROW)
  271.   110 FORMAT ( /// 4H   Z)
  272.   120 FORMAT ( /// 7H   YROW)
  273.   130 FORMAT ( ///// 41H     STEP TWO   UPDATING XROW, YROW AND Z ///)
  274.   140 FORMAT ( /// 6H   RHO // 1H , 2E13.6)
  275.   150 FORMAT (4H   R)
  276.       END
  277.       SUBROUTINE SGETRX(X,LDX,N,P,S)
  278.       INTEGER N,P,LDX
  279.       REAL X(LDX,1),S(1)
  280.       INTEGER PD2,PD2D1,I
  281.       PD2 = MAX0(P/2,1)
  282.       PD2D1 =